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Prospects for Detecting Dark Matter Halo Substructure with Pulsar Timing 
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One of the open questions of modern cosmology is the nature and properties of the dark 
matter halo and its substructures. In this work we study the gravitational effect of dark matter 
substructures on pulsar timing observations. Since millisecond pulsars are stable and accurate 
emitters, they have been proposed as plausible astrophysical tools to probe the gravitational effects 
of dark matter structures. We study this effect on pulsar timing through Shapiro time delay (or 
integrated Sachs-Wolfe (ISW) effect) and Doppler effects statistically, showing that the latter 
dominates the signal. For this task, we relate the power spectrum of pulsar frequency change to the 
matter power spectrum on small scales, which we compute using the stable clustering hypothesis, 
as well as other models of nonlinear structure formation. We compare this power spectrum with the 
reach of current and future observations of pulsar timing designed for gravitational wave detection. 
Our results show that while current observations are unable to detect these signals, the sensitivity 
of the upcoming square kilometer array is only a factor of few weaker than our optimistic predictions. 

PACS numbers: 04.50. +h, 95.36.+X, 98.80.-k 



I. INTRODUCTION 

One of the greatest puzzles of modern cosmology is the 
nature of the dark matter (DM). The latest cosmological 
observations indicate that DM has a mean cosmic mass 
density ^5 times larger than the density of the baryonic 
matter (e.g. see WMAP- 7 yr results and its presence 
is confirmed by a large amount of astrophysical evidence, 
such as rotation curves of galaxies, gravitational lensing 
effects, growth of large scale structure of the Universe, 
big bang nucleosynthesis and the dynamics of the Uni- 
verse as a whole [2j. The cosmological and astrophysical 
observations that may lead us to better understand this 
unknown component of the Universe are areas of intense 
study, focusing on DM's nature as particles, its struc- 
ture, distribution and effect on the other components of 
the Universe. Studying the small scale structure of DM, 
for example, will tell us something about the DM parti- 
cle's fundamental properties. Consequently, finding new 
footprints of DM in astrophysical observations is impor- 
tant for opening new horizons in DM studies. 

The theory of structure formation, which is based 
on gravitational instability of primordial matter density 
fluctuations and the hierarchal scheme of structure for- 
mation, assumes that collisionless DM is the main in- 
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gredient in today's cosmological structures. One of the 
features of this theory is that DM collapses into bound 
states, known as DM halos. A cold dark matter pri- 
mordial power spectrum predicts a large range of mass 
scales for these DM halos, from I0 12 — 10 14 Af© down to 
I0~ 12 — 10~ 4 AfQ Larger halos form from merger of 
smaller halos which may partly survive as substructure 
of bigger halos. 

The statistics of DM distribution and the dynamics of 
this substructure may have an effect on astrophysical ob- 
servations. One of the promising astrophysical probes for 
studying the distribution of interstellar medium (ISM) 
which is considered to be mostly baryonic matter are pul- 
sars 01 • ISM causes a dispersion on pulsars' light which 
in turn has an effect on pulsar timing residuals. 

In the present work, we push this one step further 
and study the gravitational effects of dark matter halo 
substructure on pulsar timing. This effect manifests it- 
self through the (I) Shapiro time-delay effect and 
(2) Doppler effect. The Shapiro time delay is caused by 
the presence of dark matter's dynamical potential along 
the line of sight. On the other hand, the Doppler effect 
is caused by the acceleration of the observer/pulsar be- 
cause of the pull of DM subhalos. Probing dark matter 
substructure by Shapiro time delay in pulsar timing was 
first proposed by Siegel et al. by considering the effect 
of one DM subhalo crossing the line of sight Dark 
matter studies with pulsar timing continued by Seto and 
Cooray , Pshirkov et al. [f| and recently by Ishiyama et 
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al. [^] . On the other hand the Shapiro delay effect was 
studied for relativistic neutrino and photons of SN1987A 
1C | and also for low frequency pulsars in globular cluster 
It is worth mentioning that pulsar timing was also 
used to study other astronomical effects [l2| , as a recent 
example the effect of DM subhalos crossing the line of 
sight was studied in astrometric microlensing 13 1 . 

In this work, as a complementary and more realistic 
view, we consider the statistical distribution of DM sub- 
structure and its effect on pulsars' timing residual. In 
order to study the effect of the DM substructure dis- 
tribution on pulsar timing, we need a structure forma- 
tion model. On very small scales deep into the nonlinear 
regime of structure formation, which is unaffected by halo 
merging or tidal disruption, we can use the stable clus- 
tering hypothesis. The stable clustering hypothesis was 
first introduced by Davis and Peebles [14[ as an analytic 
technique to study the galaxy correlation function in the 
deeply nonlinear regime, and was subsequently applied to 
fitting formulas for nonlinear correlation functions/power 
spectra 15, 1(|. In the current work, we use the phase- 
space stable clustering model which was recently devel- 
oped by Afshordi et al. fl7j . 

The article is structured as follows. In Sec. (II), we 
first introduce millisecond pulsars. Then in the following 
subsections we derive the power spectrum of frequency 
change of pulsars for Shapiro time delay and Dopplcr 
effects. In Sec. (Ill), we review the stable clustering hy- 
pothesis in phase-space. In Sec. (IV), we find the fre- 
quency change power spectrum and show its dependence 
on free parameters of the model, both for Shapiro and 
Doppler effects. In Sec. (V), we discuss the observational 
prospects of detecting these effects with current and fu- 
ture pulsar timing arrays. Finally, Sec. (VI) concludes 
the paper. 

For reference, we set cosmological parameters to be 
Q,° n = 0.27, er 8 = 0.8 and H a = WOh km/s/Mpc where 
h = 0.7. 



II. GRAVITATIONAL EFFECT ON POWER 
SPECTRUM OF PULSAR TIMING 

In this section, we first introduce millisecond pulsars 
as promising astrophysical observational probes to detect 
the gravitational effects of DM substructures. Then we 
derive the statistics of frequency change due to Shapiro 
and Dopplcr effects. 



A. Millisecond Pulsars 

The most stable, consistent astrophysical emitters in 
the known universe are millisecond pulsars, many of them 
remaining stable without flux change over timescales ex- 
ceeding 30 years [l8| ■ On account of this they have been 



used as precise tools to probe changes in the matter dis- 
tribution between the pulsar and earth [l9| . The pulsars 
with the highest rotational frequencies, and hence the 
shortest pulse to pulse periods, are the most stable with 
a time period of 0(1 ms). The typical residual of these 
pulsars is of order of 0(1 fj,s). This means that fluctua- 
tions in pulsar period within a short time scale (e.g. ~ 
1 hr) are less than ~ [is. These residuals do not accu- 
mulate, which means that the period remains constant 
during the time that a pulsar is stable. This is used 
to measure the pulsar's timing residuals with high accu- 
racy during a long period (~ 10 years), and to search for 
nonintrinsic changes in pulsar timing. Consequently, to 
detect any physics besides the pulsars' intrinsic changes, 
we should search for a time delay larger than the in- 
trinsic uncertainties. An important point is that many 
interesting nonintrinsic effects on pulsar timing will be 
correlated. An example is the attempt to detect gravi- 
tational waves through cross-correlation of pulsar timing 
arrays [jo| . Another possible nonintrinsic effect which 
we consider in this work is the change of the gravita- 
tional potential. The transit of DM halo substructure 
across the line of sight, which causes the Shapiro delay, 
is studied in the following subsection. This discussion is 
followed by a consideration of the Dopplcr effect, caused 
by the acceleration of pulsar/observer due to presence of 
DM substructure. 



B. Shapiro time delay 

The Shapiro time delay is caused by the presence of a 
time dependent gravitational potential along the line of 
sight. To quantify this effect, we can write the metric of 
perturbed space time as 



ds 2 



(1 + 2<S>)dt 2 + (1 - 2§)dx 2 



(D 



where <J> represents the Newtonian potential in the weak 
field limit and we set the speed of light c = 1 . In the case 
of pulsars, we can write the null geodesies for a pulse 
received at time t using the above metric as 



t = t + St 



(1 - 2<f>)dx, 



(2) 



where St is obtained from integration over the perturbed 
potential along the line of sight. As it is not possible to 
measure the absolute light travel time of any astrophys- 
ical object, we need to observe the time arrival changes 
over a detection period. The time derivative of the pulsar 
time residual is defined as 



St=- 



Sv 



(3) 



where v is the frequency of the pulsar, and Sv is the 
change in frequency (we note that this is identical to the 
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cosmological Integrated Sachs- Wolfe (ISW) effect |2lj). 
In order to find a statistical description of the induced 
time delays, we assume that DM substructures move with 
a constant velocity v across the line of sight in the in- 
direction. We refer to this as the moving screen approxi- 
mation, which enables us to relate the time derivatives to 
the spatial gradients of metric, i.e. $ = vJ^Q. The tem- 
poral correlation of the frequency changes is then given 
by 



d d 

dz^z^ 2 ( — $! — $„), (4) 



where zq is the position of the pulsar in the z-direction 
(line of sight), which we take to be ~ 1 kpc, and and 
$ n correspond to potentials at two different times. Now 
we can express the right hand side of Eq. Q in Fourier 
space as: 



4i/ 



dzidzu 



(2^)3 



d A k 



(iAf)(*^)<*(JO*(*u)>e- 



(2tt) 3 



S(5) 



By integrating over z 1 and z n , and using the definition of 
the potential power-spectrum, 



Eq. ((5]) becomes 



((-M-)n) 



d 3 ^ 



fef Zosinc 



(2tt) 3 



fcfzo 



(6) 



(7) 



where sinc(a;) = "We can take the integral over fc z , 

which results in 



<(-).(-)u> = W 



~t~tpm (8) 



x (fc x ) 2 e- ifc ^ (tl - t2 ), 

where we omit the subscript i and replace x l = utj in 
the moving screen approximation (we also assume k z ~ 
z^ 1 <C k x ,k y ). Now by considering the definition of the 
time-delay power-spectrum, 



((-),(-w4/^M^. 

V V lit J " 



(9) 



we can integrate Eq. ((5J) over the time difference of two 
observations to obtain the power-spectrum, 



^(«-)=/«-W-W« 

) v v 



where At = — £2- Now the integration over dk x and 
d (At) gives us the relation between k x and the frequency 
as k x v = ui. Consequently Eq. (|10p simplifies to: 



. . izo f dk y 



v I 2ir 




(11) 



Using the Poisson equation, we can relate the potential 
power-spectrum to the matter power-spectrum P p {k): 



(12) 



Finally by inserting Eq. (|T2l) in Eq. (|TT|) . we find the 
dimensionless uP(lu) in terms of the matter power- 
spectrum: 

uPs±(oj)\ Shapiio = 



Az Q f dk y 3 UttGY 



where we replace P p {k) = p 2 PNL,{k) , in which p is the 
mean cosmic DM density. In Sec. (IV), we will derive this 
function by using the stable clustering hypothesis. 



C. Doppler effect 

Changing the potential of DM substructure near pul- 
sars or the Earth will introduce a velocity shift, which 
affects pulsar frequencies via the Doppler effect. In the 
Doppler effect, the frequency change of a pulsar is re- 
lated to the line of sight velocity as ^ = vi. s .. So the 
correlation of the frequency changes observed at two sep- 
arate times ti and fa caused by the Doppler effect can be 
written as 



dt / dt^v.fcxV^uXu) 



<(-),(-),) = 

v v 

X e £ (*-*l) e E (*'-*2) 



where e is a small parameter to regulate the infrared 
divergence of the integral. Once more, we can write the 
right hand side of Eq. (|14l) in Fourier space. Integration 
over time variables with the limit of e — > results in 



((-M-)n> 
V V 



' k Mk)p^e- ik "< At \ (15) 



(27T) 3 ^ '{k x vf 



where we used the moving screen approximation to re- 
place time integrals by integrals over x. Again, using 
Eq. ((9|), we can write the power-spectrum of frequency 
changes as 



d(At) 



(10) 



ik x v(At) iuiAt 



I {{-U-) u )e^d{At) 
v v 

d k PJk) ^ e-i^vjAt) iuAt 



(16) 
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Integration over dk x and dAt gives us the relation be- 
tween k x and the frequency as oj — k x v. Consequently 
Eq. (fTo) results in 

If dk,, f dk, kl , . 

Because of the symmetry between the integration over 
k y and k z , we can replace k 2 y by (ky + kg)/ 2 in Eq. ([17]). 
which leads to 

p ^=lJ dk S p *{^^)> (i8) 



where k„ = \Jk'y + k^. By using the Poisson equation, we 
can relate the potential power spectrum to matter power 
spectrum, and finally write the dimensionless power spec- 
trum as 

^M| Dopplet = 

As in the case of Shapiro delay, by knowing the mat- 
ter power spectrum we can determine the pulsar timing 
power spectrum. 
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FIG. 1: Surfaces of constant average CDM phase space den- 
sity! {/)p = M?s, around a typical particle in the stable clus- 
tering hypothesis. The surfaces are assumed to be concentric 
ellipsoids (Eq. I24[) . The mass and Hubble scales at the col- 
lapse of the structure, M(£„) and H(£ a ), are related to the 
phase space density, fi£ a , on each surface via the spherical 
collapse results Eqs. (|25H28[) . while fi ~ 3% is an empiri- 
cal factor that quantifies tidal stripping and is fixed through 
comparison with numerical simulations [l7l ] . 



III. STABLE CLUSTERING HYPOTHESIS 

In this section, we discuss the stable clustering hypoth- 
esis as a model to describe the nonlinear structure for- 
mation that will give us the required matter power spec- 
trum. We make use of the phase-space stable clustering 
model recently developed by Afshordi et al. 



171. The 



collisionless Boltzmann equation at the phase-space co- 
ordinates, r + Ar, v + Av is approximately given by 



— (r + Ar, v + Av,t) 
dt 



(20) 



% + ^- ■ (v + Av) - y~ ■ [V$ + (Ar • V)V$] = 0, 
ot or ov 

where $ is the gravitational potential, and for simplicity 
we omit the vector signs of distances and velocities. We 
can reexpress the above equation in terms of the phase- 
space density in the comoving coordinates with particle 
i: 



/.(Ar, Av) = f{n + Ar, Vl + Av). 



(21) 



Using this new function, we can write the Boltzman 
Eq. (HO]) as 



df = dfi 
dt dt 



Fj .Av- J$--(Ar-V)V$ = 0. (22) 



Notice that Eq. ([22]) can be understood as the tidal limit 
of the Boltzmann equation in terms of the phase coordi- 
nates (Ar, Av), i.e. in the coordinate system comoving 
with particle i. 

The stable clustering hypothesis assumes that 
T^Hat-.Au averaged over the particles vanishes for small 
Ar and Aw. This implies that the number of neighbors 
within a fixed physical separation of a DM particle in 
the phase space does not vary with time. Now, if we as- 
sume that (/jVV$) ps (/i)p(VV$) p , then a solution to 
Eq. (jUJ) is: 



(23) 



(f) P = ^ = F ^ Ay2 + Xr - Xr! 'W' 



where F is the general solution with isotropic velocity 
distribution and N is the number of particles in the phase 
space volume of interest. By using the approximation of 
a spherically symmetric potential, the above solution can 
be rewritten by applying the Poisson equation: 



(Dp = tf. = F[(Av) 2 + 100ff 2 (6)(Ar) 2 



(24) 



<9Ar 



dAv 



where £ s and i?(£ s ) are the phase-space density and Hub- 
ble constant at the formation time of DM substructure, 
respectively (see Fig. [T]). We also use the spherical 
collapse model prediction for the halo density, which is 
roughly 200 times the critical density at the formation 
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time [22|. \i ~ 3% is the mean fraction of bound •par- 
ticle pairs that can survive the tidal disruption period, 
and is calibrated by comparison with N-body simulations 
171 ]. To determine the function F, we use the spherical 
collapse model results. The phase-space density can be 
expressed as 



ioff(6) 



(25) 



using the fact that the radius and velocity dispersion of 
halos are related as 



10Hr v 



(26) 



The phase-space volume of the collapsed halo, i.e. the 
volume of the constant-^ ellipsoid in Eq. (|23|) . is M/£ s , 
and by using Eq. (|2"5"|) , we find 



10F(6) 



10ff(6) 



(27) 



Furthermore, the mass scale that collapses at a given 
cosmological epoch is characterized by 



H 



-2/3 



<t[M(&)] ~ <5 C ~ 1.7, 



(28) 



where <5 C is the linear density threshold for the spherical 
collapse, <r[M] is the rms top- hat linear over density at 
the mass scale M, and Ho is the Hubble constant in the 
present epoch. Using the above result, the phase-space 
correlation function is obtained as 

(f(ri,v 1 )f(r 2 ,v 2 )) 

~ — / d 3 rd 3 vf(r,v)f(r + Ar,v + Av) 
Ve Jv 6 

i 

- (f(n,v 1 ))(f(r 2 ,v 2 )) + f i(f(f,v))U&r,Av). 

(29) 

In the equation above we used the assumption of ergod- 
icity to replace the ensemble average () by the volume 
average, in a given volume of phase-space Vg, while (f, v) 
are the mean values of (r"i,i>i) and {r 2l v 2 ). The second 
term is based on the stable clustering described above, 
with the assumption that \Av\ = \v\ — v 2 \ <C Avud and 
I Ar| = |n - r 2 \ < Ar tld where Aw tid and Ar tid charac- 
terize the tidal truncation radii in the phase-space. On 
the other hand, the first term in Eq. (f29|) dominates for 
large separations in the phase-space, where particles are 
not correlated. So Eq. (|29|) is an interpolation between 
the stable clustering and the smooth halo regimes. This 
is a crucial point in calculating the nonlinear power spec- 
trum of structures on small scales where it is related 
to phase-space density correlation fi(f(f,v))t; s (Ar,Av) 
term. 
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FIG. 2: Dimensionless power spectrum of density fluctua- 



k*P NL (k) 
2^ 



as a function of wavenumber k for 



tions A 2 (fc) 

the linear regime (short-dashed line), Peacock and Dodds fit- 
ting formula (dash-dotted line), Smith et al. fitting formula 
(thick dotted line), halo model (long-dashed line) and for 
k < 10 2 stable clustering hypothesis used in this work (solid 
line). 



IV. PULSAR RESIDUAL POWER SPECTRUM 
FROM STABLE CLUSTERING HYPOTHESIS 

In order to calculate the dimensionless power spectrum 
of pulsar frequency change, we need to know the power 
spectrum of matter on small scales. We now make use of 
the stable clustering hypothesis prediction, as developed 
in the previous section. 

To use the stable clustering formula obtained in Eq. 
(|29p . we must relate the matter density power spectrum 
in Eqs. (|13I19|) to the real space correlation function of 
densities. In the stable clustering hypothesis in phase- 
space, on small scales this relation becomes 

(pftWu)) = f &%&<ru(f (so) 

(filcfAvfiifir, 0)>&(Ar, Av) 

= UPavg J d 3 Av£ s (Ar,Av). 

In order to find the dependence of the dimensionless 
power spectrum ujPsv_(ui) on uj, we should have the rms 
top-hat linear overdensity a(M). a(M) is the integral of 
linear matter power spectrum on a chosen window func- 
tion as 



a\M) 



d 3 k 
(27r) a 



P L (k)W 2 {kR), 



(31) 
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where -Pl(^) and W(kR) are the linear matter power 
spectrum and the Fourier transform of the spherical top- 
hat filter of radius R, respectively, where 



P L (k) = Ak n °T 2 {k), 



W(x) 



3(sinx — xcosx) 



(32) 



(33) 



Here n s is the scalar spectral index of primordial mat- 
ter power spectrum and M = 47ri? 3 p m /3. The transfer 
function can be approximated by the BBKS [23[ fitting 
formula, 



T(k = qfl m h z Mpc- L ) 



2* , -is „ Ml + 2-34g] 



2.34<j 



(34) 



[1 + 3.89q + (16.2q) 2 + (5A7qy + (6.71 9 ) 4 ] 



41-1/4 



Using Eqs. (|30I31|) . and the stable clustering hypothe- 
sis, we can find an expression for the matter power spec- 
trum on small scales, shown in Fig. ^ for a mass range 



10 6 to 10 12 solar masses and ^ = 0.03. For qualitative 
comparison, we also plot the dimensionless power spec- 
trum A 2 (k) = l (k) obtained from the halo model 
of structure formation [24f , as well as the fitting formu- 

16] and Smith et al. for the 



las of Peacock and Dodds 
nonlinear power spectrum 



25} . We note that these ap- 



proximations are based on fits to numerical simulations 
at k < 10 2 Mpc -1 , while the stable clustering hypothesis 
(which goes into second term in Eg. ([2"9"|) 'l is expected to 
hold for k ^ 10 Mpc -1 , and thus should be a more appro- 
priate measure of small scale dark matter structures. On 
larger scales, the matter power spectrum is dominated 
by the first term on Eq. (|2"5|) . which is equivalent to the 
standard halo model (i.e. the long-dashed line in Fig. @). 
The cut-off in the halo model power spectrum is related 



to the size of smallest halo mass of M„ 



10~ b M, 



0- 



Now, using the nonlinear power spectrum obtained 
from stable clustering, the dimensionless power spectrum 
for the Shapiro time- delay effect can be written as 



n / \ i 4z o f dk v 3 (4ttG) 2 f A .. .sin(fcAr) f A ,„ , A . , 10H[£ s (Ar, Av)] . . 

In the case of the Doppler effect, we can also find the dimensionless power spectrum of pulsar frequency change in 
terms of phase-space density derived from the stable clustering hypothesis. In this case Eq. (|19[) is expressed as 

UPhaio f dk v k z y (4ttG) 2 /"^A^^A^sin^Ar) f 2 10/f[£ s (Ar, Av)] 



WP *^I- = = / ^-k^ X J MAr)'d(Ar) n ^ J d(Av)4ir(Av) ^2 Af ^ 8 ( Ar , Av)] ' (36) 



Notice that k = J (^) 2 + fc 2 , and the Hubble parameter 

and the mass are related by a(M) through Eq. (|28p . phaio 
is the smoothed halo local density at solar system which 
is assumed to be ~ 10 5 p cr it- 

In order to numerically perform the third integrations 
in Eqs. (|35H36[) over (k y ,Ar,Av), we trade Av with 
M as the integration variable using F^ 1 = (Av) 2 + 
100i? 2 (Ar) 2 (see Fig. [T]). We then perform the integra- 
tion in three steps: 

1. Noting that fixing M and Ar fixes Av through 
Eqs. ()24II28|) . we can first perform the k y integral 
for fixed Ar and Av. Since the integrand can have 
fast oscillations in k y , we find asymptotic expan- 
sions for the k y integral in the wAr/v ^> and <C 1 
limit, and devise an interpolation between the two 
regimes with less than 1% error, compared to the 
exact integral. 

2. We then perform the Ar integral from zero to the 
maximum of ( J F 1 - 1 ) 1 / 2 /(10i7), which is fixed by 



mass M, through Eq. (|2"4"]) . 

3. Finally, we take the integral over substructure 
mass, M, from M m i n to M max , which we discuss 
below. We note that H[£ s (Ar, Av)] also becomes a 
function of M. 

Now we are able to calculate the dimensionless power 
spectrum wPfo (w) in terms of the frequency u) numeri- 
cally, for Shapiro time delay and the Doppler effect. For 
convenience we define the dimensionless parameter h p as 

h p = [^-u,P(u,)]i, (37) 

which is shown in Fig. ^ for the Doppler effect (solid 
line), Shapiro delay (dot-dashed line), and white noise 
(dashed line, which is computed in Appendix [A"|) . In 
order to calculate h p numerically, we consider a realistic 
set of parameters (but later study the effect of changing 
these parameters) . We choose the velocity of dark matter 
substructures v = 300 km/s (typical of relative velocities 



7 



Timescale (year) 
10' 10' 





i 1 1 1 1 


ill l 


I | I I I I I I I I 


- 




Doppler 




: 
- 




Shapiro 









White noise 






V= 300 km/s 


H=0.03 


1|is 


- 

: 
- 


M = 10" 6 M 


z=lkpc 


^ ■ — ' 1 00 ns^ 


_ 

: 
■ 

- 




^ ■ — * 









































































































~i 


1 i i i i 


1 i i i 


I i i i i I i i i i I 



-9.5 -9 -8.5 -8 -7.5 



log#/Hz) 

FIG. 3: Pulsar residual power spectrum as a function of fre- 
quency (bottom x-axis) and the span of observation time in 
years (top x-axis) for time delay caused by the Doppler effect 
(solid line) and time delay caused by Shapiro effect (dash-dot 
line). The long dashed lines represent levels of white noise for 
100 ns (bottom) and 1 [ms (top) measured biweekly (see the 
Appendix) 

in the Milky Way halo), the typical distance of pulsars 
to Zo = 1 kpc, the mean fraction of bound particles that 
can survive the tidal disruption period fx = 0.03 [l7j], the 
minimum mass of DM substructure M min = 1O~ 6 M 
and also the maximum M max — 10 12 M Q (the total mass 
of a galactic halo). Later we will show that h p is almost 
independent of M rnax . 

The power spectra of Shapiro and Doppler effects in 
Fig.© are well described by power- laws: 

^MI Shapiro tx wP^H| D<wtar cx (38) 

These behaviors can be understood by noticing that the 
Av integral (i.e. the last integral) in Eqs. (|35H36j) scales 
as i? 2 (£ s ), if we use the spherical collapse relations of 
Sec. piip . Since most small structures with CDM initial 
conditions collapse around the same time, this is approx- 
imately constant. The contribution to the rest of the 
integrals is dominated by ky 1 ~ Ar ~ v/lu, so the inte- 
gral over distances scales as (Ar) 3 cx w~ 3 . Plugging this 
into Eqs. (|35H36[) yields the scalings of Eq. (|38| . 

To physically understand the scaling for the Doppler 
effect we can once more Fourier transform the power- 
spectrum in Eq. (f36)) to find that voop. ~ ^7 is propor- 
tional to vt 2 = (vt) x t, i.e. the magnitude of acceleration 
is proportional to distance traveled by the earth/pulsar. 
This is exactly what one expects for the gravitational 
field in a medium with roughly uniform density, and is 
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10' 10° 
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FIG. 4: Pulsar residual power spectrum as a function of fre- 
quency (bottom-x axis) and the span of time in years (top 
x-axis) for time delay caused by the Doppler effect (top-line) 
and time delay caused by Shapiro effect (bottom line) for a 
maximum mass of halo M max — 10 12 Mq (dash-dot line) and 
M max = 1O 8 M (solid line). 

due to the fact that most small substructure forms at 
roughly the same density cx -ff 2 (£ s ). However, the direc- 
tion of acceleration is random, as different substructures 
will dominate the local gravity on different scales. 

An important point to consider before examining the 
effect of different parameters on pulsar timing is the 
study of the effect of maximum mass in the integrals. 
As we show in Fig. (j4]), the total dependence of h p 
on maximum mass is small, where we plot the h p for 
Mmax = 10 12 M Q , the total mass of a typical galaxy and 
M max = 10 8 A/ Q , for a more realistic tidal cut-off for sub- 
haloes at our position in the Milky Way. This confirms 
that, not surprisingly, most of the observable effects on 
pulsar timing comes from CDM small scale structure. 

Now we examine the dependence of the power spec- 
trum on different parameters of the model. We plot the 
dimcnsionlcss amplitude h p for the Doppler effect for dif- 
ferent velocities of dark matter substructures and the fi- 
parameter of stable clustering in Fig. ([5]), which shows 
that h p is proportional to velocity and the square root of 
the 11 parameter. 

In Fig. ([6]), we plot the power spectrum for different 
mass minima of DM substructures. As shown in Fig. ([6]), 
the uj~ 4 dependence of h p does not change by changing 
the minimum of the mass. However, the amplitude of 
the signal increases when the interval of integration is 
increased. 

In Figs. and © we plot h p given different primor- 
dial spectral index n s , for Doppler and Shapiro effects 
respectively. For n s < 1, the slope of h p does not change, 
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FIG. 5: Pulsar residual power spectrum as a function of fre- 
quency (bottom-x axis) and the span of time in years (top 
x-axis) for Doppler effect (solid line) for different velocities 
and fi of dark matter substructures. 



FIG. 7: Pulsar residual power spectrum as a function of fre- 
quency (bottom-x axis) and the span of time in years (top x- 
axis) for Doppler effect for different primodial index of matter 
power spectrum. 




FIG. 6: Pulsar residual power spectrum as a function of fre- 
quency (bottom-x axis) and the span of time in years (top 
x-axis) for Doppler effect (solid line) for different minimum 
masses of dark matter substructures. 



FIG. 8: Pulsar residual power spectrum as a function of fre- 
quency (bottom-x axis) and the span of time in years (top x- 
axis) for Shapiro effect for different primordial index of matter 
power spectrum. 



as a(M) becomes flat for low masses. On the other hand 
for larger n s , we see a shallower u> dependence for h p as 
there is more power on small scales. 

V. OBSERVATIONAL PROSPECTS 

Finally, to explore the observational prospects for the 



detection of pulsar frequency change due to dark matter 
substructures, we compare our results with the obser- 
vational bounds put on detection of gravitational waves 
(GW) by pulsars. The observed quantities are similar in 
both cases, and the power spectrum of pulsar frequency 
change caused by Doppler or Shapiro effects is red sim- 
ilar to gravitational waves. That is, there is an excess 
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power at low frequencies, or long timescale correlations 
in residuals. The gravitational wave effect on pulsar tim- 
ing is also in the nHz frequency range 2(J 26| . similar 
to the substructure effect we are considering here. 

In particular, the frequency change due to gravita- 
tional waves is roughly ~ hy, the amplitude of gravi- 
tational waves, which allows us to directly translate con- 
straints on hij, to constraints on 8v/v. Moreover, similar 
to the characteristic quadrupolar pattern that gravita- 
tional waves induce in pulsar timing residuals {e.g., [Ifij]), 
the Doppler effect induces a dipolar pattern in the sky, 
which can be used to distinguish it from intrinsic changes 
in individual pulsars. 

More specifically, the frequency shift due to the com- 
bination of Doppler effect, gravitational waves, and in- 
trinsic effects is given by: 

a = I a (t) + n a • v(t) H F , (39) 

v 1 + k • n a 

where n a is the unit vector along the direction of pulsar 
a, I a {t) is the frequency shift intrinsic to the pulsar, v(t) 
is the earth's velocity, and h(k, t) is the amplitude of a 
gravitational wave with wave-vector k. The cross-power 
spectrum of frequency-change between different pulsars 
is given by 

Psz.(w)\ab = P^(u)\int.S(Xab) 
+ Psu (uj)\ Doppler (1 - 2x ab ) + P(u>) | grav . c{x ah ) ,(40) 



where 



and 



x ab = (1 - n a ■ h b )/2, 



c(x) 



-a; In a; 



x 
1 



(41) 



(42) 



is the expected correlation pattern of timing residuals 
for an isotropic stochastic gravitational wave background 
[26| . Therefore, pulsar timing cross-power spectra are af- 
fected by the intrinsic, Doppler and gravitational waves, 
P^{w)\int.,P$v(w)\Doppler, and P(u)\ g rav. with different 
angular dependences, which can be used to distinguish 
these effects. 

In Fig. (|9|) we plot the realistic and optimistic predic- 
tions for detection of h p , which is similar to the gravi- 
tational wave dimcnsionlcss strain, and compare it with 
the current observational limits from a pulsar timing ar- 
ray [3] considering the sensitivity limit for time residuals 
of observed millisecond pulsars obtained via (see the ap- 
pendix of [27J for details) 



h l : m 



x 



fitrmsf 



iV p 1/2 (TA/) 1 /4 : 



(43) 



where hp™ 1 is the sensitivity limit of detectors, 6t rms = 
yj (St 2 ) is the root mean square value of the timing resid- 
uals, A/ is the frequency bandwidth of search, N p is 
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FIG. 9: Pulsar residual power spectrum as a function of the 
frequency for Doppler effect for realistic and optimistic sig- 
nals (see text for definition of realistic and optimistic param- 
eters). The limits from current and future experiments are 
also shown. 



the number of pulsars, and T is the time span of ob- 
servation. The pulsar timing array sensitivity is scaled 
with frequency as hp 1 ™ 1 cx / and reaches a minimum at 
a detectable frequency of / ~ 1JT. This produces the 
wedge-like sensitivity limit curves in Fig. @. The sensi- 
tivity limit is also proportional to St rms , improving as the 
precision of pulsar timing residuals detection is increased. 
By increasing the observational time, we increase the sen- 
sitivity and also the span of frequency. 

We also plot the predicted sensitivity of Parkers pul- 
sar timing array (PPTA) [28[ and the square kilome- 
ter array (SKA) [29( for h p . The upper bounds for 
future PPTA and SKA experiments are obtained from 
the detectable time residual correlation of simulated pul- 
sars with consideration of all instrumental, calibration 
and observational errors (such as pulsar intrinsic period 
changes and glitches) • For example the PPTA bound 
is obtained by considering 20 radio pulsars for 5 year 
with 5t rms = IOOtis which provides a peak sensitivity of 
h lim ~ 2 x 10~ 15 at / w 7 x 10~ 9 . For SKA, with the 
same number of pulsars, the sensitivity is improved by 
increasing the span of observation to 10 years with tim- 
ing accuracy 5t rms = 10ns, leading to a constraint on 



the pulsar residual power spectrum of 
/« 7 x 10~ 9 



1.6 x 10 



-16 



at 



27|. 



Finally, we study the effect of uncertainty in the models 
of nonlinear structure formation on our results. In other 
words, how much will our results depend on the choice of 
stable clustering hypothesis? As we argued above, sta- 
ble clustering is the only known physical prediction for 
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FIG. 10: Pulsar residual power spectrum as a function of 
the frequency due to Doppler effect, for realistic signals in 
different models of nonlinear structure formation. 



the nonlinear power spectrum on very small scales. Nev- 
ertheless, we can calculate h p for the nonlinear power 
spectra of other clustering models in Fig. ©. The rela- 
tive magnitude of h p in two different models is obtained 
from Eq. ([T9l) : 



(ml) 

Lp 

Am2) 



au _Lp( mi : 

k 4 



(*))/( / dk 



K p (m2) 



NL 



(*0) 



1/2 



(44) 

where superscript (ml) and (to2) indicate the models. In 
Fig. (|TU|) . wc plot h p for different models of nonlinear 
structure formation by using the realistic parameters for 
the models. 

An interesting point to notice is that different mod- 
els of nonlinear structure formation have (almost) the 
same frequency dependence, LoP^ v j v oc a; -4 , as in stable 
clustering. This is because of the moving screen approx- 
imation k x v = uj and k z ~ z _1 -C k x ,k y , which is ap- 
plicable in the alternative models as well. On the other 
hand, the main contribution of the integrals in Eq. (|4"4")l 
from the nonlinear matter power spectrum comes when 
k y ~ lo/v ~ 10" 8 Hz/300 km - 10 9 Mpc" 1 . In this case 
h™ 1 /h™ 2 reduces to the ratio of Pnl's, which is nearly 
independent of wavenumber (and thus frequency; see Fig. 
[2]) for relevant scales . 

Closer examination indicates that the frequency de- 
pendence of the Smith et al. model is slightly shallower 
than the others (hp mlth ~ w -1 - 93 ). This is due to the 
fact that Smith et al. predict a much bluer spectrum on 
small scales ( Figj2]), which is similar to the case of stable 
clustering with higher power index (FigJT])- 

In summary, we find that the signature of the Doppler 



effect in pulsar timing is largely independent of the non- 
linear structure formation model, which only introduces 
(a factor of a few) uncertainty in the amplitude of timing 
residuals, h p . Our results show that while current ob- 
servations are unable to detect the effect of dark matter 
substructure on pulsar timing, error projections for the 
upcoming square kilometer array (SKA) are only a factor 
of few higher than our optimistic predictions. 

In the end, we should note that, unlike the Doppler 
effect, the Shapiro time delay does not have a coher- 
ent pattern on the sky, as different lines of sight are 
largely uncorrelated. This makes it much harder to dis- 
tinguish Shapiro time delay from pulsar intrinsic fre- 
quency changes. 

VI. CONCLUSIONS AND DISCUSSION 

In this work, we studied the gravitational effect of 
DM substructures on pulsar timing, through Doppler 
and Shapiro (or ISW) effects. We calculated the dimen- 
sionlcss power-spectrum of a pulsar's frequency-change, 
which is related to the matter density power-spectrum in 
the nonlinear regime. We used the stable clustering hy- 
pothesis to extract the nonlinear matter power-spectrum, 
and showed that the frequency-change is dominated by 
the Doppler effect. Next we varied the free parameters 
of the model, which had the following effects on the di- 
mensionless power, h p : 

1. h p due to Doppler effect is linearly proportional to 
velocity of DM substructures. 

2. The main contribution of DM substructures comes 
from the minimum mass in DM hierarchy: as we 
increase the domain of integration over DM subhalo 
masses, we get more signal. 

3. hp has a dependence on fj, 1 / 2 , the fraction of particle 
pairs that remain bound in the stable clustering 
hypothesis. 

4. hp due to the Shapiro effect scales as the square 
root of distances to pulsars, as it depends on the 
integrated gravitational effect over the line of sight. 

5. For larger primordial spectral index, n s , the fre- 
quency dependence of h p is shallower, because the 
main contribution of h p comes from low masses, 
where the power is increased. However, for n s < 1, 
the frequency dependence becomes independent of 



6. The frequency dependence of h p is nearly indepen- 
dent of nonlinear structure formation model, al- 
though its amplitude could change by a factor of 
a few. 

Finally, we compared the dimensionless power spectrum 
of pulsar frequency change, for realistic and optimistic 
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sets of parameters, with current and future pulsar tim- 
ing experiments, designed for detection of gravitational 
waves. Our results show that our optimistic estimate of 
the hp signal is only a factor of a few smaller than the 
sensitivity of the planned square kilometer array (SKA), 
making this method a potentially promising avenue for 
the detection of DM substructure on very small scales. 
While this may sound too futuristic, it is worth noting 
that more dedicated pulsar timing follow-ups of pulsars 
discovered by SKA, as well as better noise removal tech- 
niques for ISM contamination of timing signals (e.g., (3l| ) 
will be able to potentially push down the noise below our 
conservative forecasts. 

We should further note that, as can be seen in Fig. 
(p~0|) . if this signal is ever detected, there will be degenera- 
cies between parameters that quantify the nature of DM 
and those of structure formation (in both linear and non- 
linear regimes). Therefore, further study into the nature 
and properties of the signal (or independent observables) 
will be necessary to disentangle these degeneracies. 

While our paper lays the groundwork for future sta- 
tistical detection of dark matter substructure through 
pulsar timing, many practical challenges and theoreti- 
cal uncertainties remain. Here we point out two, along 
with potential resolutions: 

First, it is important to note that the observed Doppler 
effect in pulsar timing depends on the total gravita- 
tional acceleration, which can be contributed by nearby 
stars/planets, in addition to local dark matter substruc- 
ture. However, the gravitational pull of stars/planets on 
Earth can be calculated by knowing their masses and po- 
sitions around Earth, and thus, in principle, can be com- 
puted and corrected for (e.g., |32|). Similar effects on 
the acceleration of pulsars will be uncorrclated for differ- 
ent pulsars, and thus can be distinguished from Earth's 
acceleration. 

A second concern is the possible non-Gaussianity of 
the signal. For example, microlcnsing events due to 
stars in the Galactic halo could lead to large magnifica- 
tions, but have very small optical depth, and thus hap- 
pen rarely. Therefore, the power spectrum gives a very 
incomplete description of the observables in microlens- 
ing events. However, in contrast to magnification events 
that trace projected density, the gravitational effects on 
pulsar timing that we discuss here trace the integrated 
potential, which is much more smooth. Moreover, the 
small CDM substructure is much more diffuse than stars, 
which further reduces the skewness of the signal. There- 
fore, unlike microlensing events, the observed signal is 
likely to be contributed by a variety of structures on dif- 
ferent scales (e.g. Fig. |6|) with no sharp boundaries. 
This is why we expect a close to Gaussian signal, simply 
based on the central limit theorem, which suggests that 
the power spectrum might provide adequate statistical 
description of these effects. 
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Appendix A: Statistics of a z for stability of Pulsars 
and White noise calculation 



The ability of a pulsar timing array to detect any de- 
lay in the received pulses to measure the dark matter 
halos substructures depends on the pulsar timing stabil- 
ity. Timing stability is related to how long the rms of 
timing residuals can be kept small, from which we can 
estimate the potential to detect Doppler and Shapiro ef- 
fects. Statistical artifacts such as a large gap in data 
sampling, or a large variation in error-bar size, may pre- 
vent a reliable power spectrum of pulsar timing data. An 
alternative approach is a z statistics, as described by e.g. 
Matsakis et al. 13311: 



„2\l/2 



(Al) 



where () denotes the average over subsets of the pulsar 
timing data, and C3 is determined from a polynomial fit 



c + ci(i - t ) +c 2 (t- t ) 2 +c 3 (t- t f 



(A2) 



to timing residuals for each subset, and r is the length 
of the subsets. In order to connect our theoretical calcu- 
lations to the observed pulsar time residuals we should 
find a relation between a z and the calculated power spec- 
trum. From the polynomial fit to the timing residuals we 
find that 



1 d . •• , 1 d bv 

c 3 ~ At s ~ , 

6 dr 6 dr v 



(A3) 



where we assume that the fitting procedure depends on 
At | s , which is coarse grained on the scale of r. The 
correlation function C3 can be written as 

(4) = ±{({-)V(i-)U-)\J}- (A4) 
18 v v v 

Now, using Eqs. (|9IAllA4|) we obtain 
a z (r)^^\J^ ^P(u)[1-cos(ut)]\ . ( A5) 
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In order to find the white noise corresponding to pulsar 
timing we derive the relation of the dimcnsionlcss power 
spectrum of pulsars with the sampling time and the un- 
certainty in the pulsar timing measurement. The cross 
correlation of time residuals of pulsar timing is related to 
the accuracy of measurement t a as 



(St(h)St{t 2 )) = {t a ) 2 S tlt2 , 



(A6) 



where <5t 1 t 2 is the Kroneckcr delta and St is the time resid- 
ual of pulsar timing related to frequency change as 



St 



5v 



-dt 



(A7) 



Now the power spectrum of time residuals is obtained as 



t (St(t 1 )St(t 2 ))dt = Tt 2 a , (A9) 



which yields the dimcnsionlcss power spectrum, 



h p = 



bjPdv (uj) 



1/2 



. ,3/2 



(A10) 



The correlation of timing residuals can be approximated 
in the time span of r, which is the period of sampling as: 



(St(h)St(t 2 )) ~ (tafrSih - t 2 ) 



(A8) 



To find the white noise lines in Fig- d3j) , we set the sam- 
pling time of pulsar timing r to be 2 weeks and the ac- 
curacy of pulsar timing, t a , to be 100ns and 1/xs. 
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